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Abstract 

Serial  simulation  is  required  to  predict  the  behavior  of  an  electrochemical  system  undergoing  many  processes.  This  is  demonstrated 
through  the  simulation  of  charge/open-circuit/discharge  processes  of  a  thin  film  nickel  hydroxide  electrode.  The  numerical  issues  involved 
in  this  kind  of  simulation  are  discussed.  An  efficient  and  robust  procedure  is  presented.  It  can  be  easily  used  to  achieve  the  serial  simulation 
of  electrochemical  processes,  e.g.,  battery  cycling,  cyclic  voltammetry  etc.  ©  2001  Elsevier  Science  B.V.  All  rights  reserved. 
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1.  Introduction 

Modeling  of  electrochemical  systems  usually  yields  large 
sets  of  differential  and  algebraic  equations  (DAEs)  or  partial 
differential  and  algebraic  equations  (PDAEs).  Since  PDAEs 
can  always  be  converted  to  DAEs  by  some  approximation 
methods,  e.g.,  the  method  of  lines  (MOL)  [1],  only  DAEs 
will  be  discussed  here.  DAEs  can  be  generally  expressed  in 
an  implicit  form: 

F(t,y,y')=  0  (1) 

where  F,y,y'  e  Rn.  A  system  of  DAEs  is  characterized  by 
its  index,  which  is  defined  as  the  minimum  number  of 
differentiations  to  convert  the  DAEs  into  equivalent  ODEs 
[2-4].  ODEs  can  be  regarded  as  special  DAEs  with  index 
equal  to  0.  DAEs  with  index  higher  than  1  are  usually 
difficult  to  solve  and  are  under  active  research  at  present. 
DAEs  with  index  equal  to  1  are  normally  encountered  in  the 
modeling  of  electrochemical  systems.  They  behave  simi¬ 
larly  to  stiff  ODEs,  and  are  generally  solved  with  similar 
implicit  methods.  One  popular  approach  to  solve  index- 1 
DAEs  is  the  backward  differentiation  formulae  (BDF) 
method,  which  has  been  implemented  in  many  DAEs  solvers 

[2,3]- 
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Simulating  a  single  process  of  a  physical  system  (i.e., 
solving  one  set  of  ODEs/DAEs)  is  straightforward.  How¬ 
ever,  an  electrochemical  system  (e.g.,  a  battery)  rarely 
operates  for  one  process  only;  instead,  it  usually  undergoes 
several  different  processes  in  series.  To  predict  the  corre¬ 
sponding  behavior,  simulation  of  many  processes  consecu¬ 
tively  is  needed,  which  has  been  researched  in  recent  years 
as  the  simulation  of  combined  continuous/discrete  processes 
[4]  and  the  simulation  of  hybrid  systems  [5].  This  kind  of 
simulation  is  demonstrated  here  for  a  nickel  hydroxide 
electrode. 

Fig.  1  presents  a  schematic  diagram  of  a  thin  film  nickel 
hydroxide  electrode.  Modeling  of  the  electrode  behavior 
yields  the  following  governing  equations: 


<9(1  -yi) 

dt 


=  Dn+ 


<92(l  -yt) 


dx2 


with  boundary  conditions 


Dn+ 


<9(1  -  yi) 

dx 


-Dn+ 


P  <9(1  -  yi)  =h 

W  dx  x=l  F 


(2) 


(3a) 

(3b) 


and 


ji  T  ji  4Pp  —  0 


(3c) 
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Nomenclature 

Dh-  proton  diffusion  coefficient  in  the  nickel  active 
material  (cm2/s) 

F  Faraday’s  constant  (96  487  C/eq.) 

F  differential  algebraic  equations 

F(k)  differential  algebraic  equations  for  Process  k 

g  equations  for  discrete  events 

g(k)  equations  for  discrete  events  for  Process  k 
i  index  for  discrete  nodes 

zapp  applied  current  density  on  the  nickel  electrode 
(A/cm2) 

z'oi  exchange  current  density  of  the  nickel  reaction 

(A/cm2) 

z02  exchange  current  density  of  the  oxygen 
reaction  (A/cm2) 

j  i  current  density  of  the  nickel  reaction  (A/cm2) 

j2  current  density  of  the  oxygen  reaction  (A/cm2) 

/  thickness  of  the  nickel  active  material  (cm) 

N  number  of  nodes  used  in  the  spatial  discretiza¬ 

tion 

Pi  interpolation  polynomial  for  dependent  vari¬ 

able  yt 

R  ideal  gas  constant  (8.3143  J/mol/K) 

Rn  n-dimensional  vector  in  real  domain 

t  independent  time  variable  (s) 

t(k)  starting  time  for  Process  k  +  1  (s) 

T  temperature  (K) 

W  molecular  weight  of  Ni(OH)2  (g/mol) 

x  spatial  coordinate  across  the  film  of  the  nickel 

active  material  (cm) 
y  dependent  variables  in  DAEs 

y'  time  derivatives  of  dependent  variables  in 

DAEs 

yx  mole  fraction  of  NiOOH 

y2  potential  difference  at  the  solid-liquid  inter¬ 

face  (V) 

Greek  letters 

</>eq ;1  equilibrium  potential  of  nickel  reaction  (V) 

^eq,2  equilibrium  potential  of  oxygen  reaction  (V) 
p  density  of  the  nickel  active  material  (g/cm3) 


where 


h  =  *'oi  |2(1  -  yi)  exp 0^  (y2  -  </>eq;1)^ 

-  2yi  exp  (y2  -  </>eq,i))  ]  (4a) 

h  =  *02  [exp  (y2  -  (J)eq  2)j  -  exp  ^  (y2  -  </>eq,2))  | 

(4b) 

The  initial  conditions  are 


yx(0<x<  l,t  =  0)  =0.05  (5a) 

y2(t  =  0)  =  0.40  V  (5b) 


In  the  above  equations,  yi  is  the  mole  fraction  of  NiOOH 
and  y2  is  the  potential  difference  at  the  solid-liquid  interface. 
Parameter  values  are  listed  in  Table  1.  The  value  of  zapp 
varies  for  different  processes,  i.e.,  it  has  positive  values, 
zero,  and  negative  values  for  charge,  open-circuit,  and 
discharge  processes,  respectively.  To  solve  the  above  model 
equations,  the  MOL  is  used  first  to  convert  PDAEs  to  DAEs, 
i.e.,  the  spatial  derivatives  in  Eqs.  (2)  and  (3)  are  approxi¬ 
mated  with  three-point  finite  difference  formulae  on  N 
uniformly  spaced  nodes: 

<9(1  -  yi  [*])  _  ^  -yi  [*  -  1]  -  yi  [*'  +  1]  +  2yi  \i] 

— rn —  ~  Dh - x? - 

for  1  <  i  <  N  (6a) 


-PH+  3yi  1*1  ~  4yi  121  -  — 131  =  0  for  i  —  0 
Ax 


(6b) 


p  — yx [A*  —  2]  +  4yi [A  —  1]  —  yi[N\  _  j\ 
n+W  Ax  F 


h  +h  ~  4PP  =  0  for  i  =  N 


for  i  —  N 
(6c) 
(6d) 


where 

^=1^  (7) 

The  DAEs  in  Eq.  (6)  have  A+l  time-dependent  variables 
where  N—2  of  them  Cy  1  [/],  1  <i<N)  are  differential  variables 


Nickel  Active  KOH 

Material  Electrolvte  Solid-Liquid 


Fig.  1.  Schematic  of  a  film  nickel  hydroxide  electrode. 


Table  1 

Parameter  values  in  the  model  of  the  nickel  hydroxide  electrode 


Parameter 

Value 

T 

298.15  K 

‘/’eq.l 

0.420  V 

(hq,2 

0.303  V 

P 

3.4  g/cm3 

W 

92.7  g/mol 

l 

lxl0~4  cm 

Du 

5xl0~12  cm2/s 

ioi 

lxlO^4  A/cm2 

4)2 

lx  10~10  A/cm2 
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and  three  of  them  (ji[l],  ydN]  and  y2)  are  algebraic 
variables. 

For  the  above  model,  the  following  processes  need  to 
be  simulated: 

•  Process  1:  charge  at  /app=l  x  10-4  A/cm2  for  3600  s  or 
to  the  cutoff  voltage  0.60  V; 

•  Process  2:  open-circuit  (7app=0)  for  1800  s  or  to  the 
cutoff  voltage  0.30  V; 

•  Process  3:  discharge  at  iapp=— lxlO-4  A/cm2  for 
3600  s  or  to  the  cutoff  voltage  0.20  V; 

•  Process  4:  open-circuit  (/app=0)  for  1800  s  or  to  the 
cutoff  voltage  0.30  V. 

In  addition  to  solving  the  governing  equations  given 
above  for  the  continuous  processes,  the  simulation  has  to 
handle  the  specified  termination  conditions  for  each  process, 
which  can  be  expressed  as  follows: 


Process  1: 

t  -  t(0)  -  3600  =  0 

(8a) 

y2  —  0.6  =  0 

(8b) 

Process  2: 

o 

II 

o 

o 

00 

1 

1 

(9a) 

y2  —  0.3  =  0 

(9b) 

Process  3: 

t  -  t(2)  -  3600  =  0 

(10a) 

y2  —  0.2  =  0 

(10b) 

Process  4: 

t  -  t(3)  -  1800  =  0 

(11a) 

y2  —  0.3  =  0 

(lib) 

The  termination  specifications  of  a  continuous  process 
have  been  called  discrete  events  [4,5].  The  equations  of  the 
discrete  events  are  independent  of  the  governing  equations 
of  continuous  processes.  In  the  above  example,  the  discrete 
events  are  given  by  the  time  specification  and  also  by  the 
cutoff  value  of  potential  y2  for  each  process.  The  occurrence 
of  a  discrete  event  terminates  the  current  process  and  starts 
the  next  process. 

2.  Serial  simulation  of  different  processes 

To  simulate  the  processes  given  above  in  series,  a  coupled 
system  of  DAEs  and  discrete  events  must  be  handled  [4],  as 
shown  in  Fig.  2  and  listed  below: 

(12a) 

II 

© 

(12b) 

Fm(t,y,y')  =  0,  [/<V2>] 

(12c) 

Fig.  2.  Schematic  of  the  simulation  of  a  series  of  processes. 


£(2)(r,3;)  =0,  [*W,*(2)]  (12d) 

FW(/.j./)=0.  [t(n-|),t(”)]  (12e) 

g(n)(t,y}^0,  (12f) 

Several  numerical  issues  need  to  be  addressed  to  solve  the 
above  system  of  equations  [4].  The  first  issue  is  the  con¬ 
sistent  initialization  of  a  continuous  process.  For  ODEs,  the 
initialization  setting  is  trivial,  i.e.,  specifying  the  values  of 
dependent  variables.  However,  for  DAEs,  the  freedom  of 
initial  settings  is  less  than  the  number  of  equations,  and 
specifying  a  consistent  initialization  is  usually  a  nontrivial 
task.  If  inconsistent  initialization  values  are  provided,  many 
DAEs  solvers  will  fail  on  the  first  integration  step  [6,7].  To 
achieve  consistent  initialization  for  DAEs,  a  widely  used  ad 
hoc  approach  is  to  take  an  implicit  Euler  integration  for  a 
very  small  time  step  without  continuity  constraint  on  the 
algebraic  variables.  With  this  approach,  the  continuity  of 
differential  variables  is  maintained,  while  algebraic  vari¬ 
ables  may  have  jumps  in  their  values,  which  is  undesirable 
when  algebraic  variables  are  more  accurately  known.  In  this 
work,  a  flexible  and  robust  initialization  solver  DAEIS  [7] 
has  been  utilized. 

DAEIS  is  designed  for  index- 1  DAEs  with  differential  and 
algebraic  variables  explicitly  identified,  which  is  encoun¬ 
tered  frequently  in  the  modeling  of  electrochemical  systems, 
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and  can  be  expressed  as 

^,ji,y1,j2)  =  o  (i3) 

where  y  i  and  y2  are  vectors  of  differential  and  algebraic 
variables,  respectively.  Supposing  the  DAEs  in  Eq.  (13) 
consist  of  n  equations,  p  differential  variables,  and  n—p 
algebraic  variables,  to  determine  uniquely  all  initial  values 
of  the  dependent  variables  and  time  derivatives,/?  degrees  of 
freedom  exist.  In  DAEIS,  the  calculation  of  a  consistent  set 
of  initializations  for  Eq.  (13)  is  treated  as  specifying  p 
variables  of  ji,  y2  and  y\  and  solving  Eq.  (13)  for  the 
remaining  n  variables  of  jq,  y2  and  y\ ,  which  is  essentially 
a  nonlinear  equation  solving  problem.  DAEIS  allows  the 
initial  values  of  accurately  known  variables  (including  alge¬ 
braic  variables)  to  remain  unchanged  and  the  initial  values 
of  other  variables  to  be  determined  during  the  initialization 
process. 

Consistent  initialization  is  not  only  needed  for  the  first 
process.  When  switching  between  two  different  processes, 
inconsistency  of  initial  values  naturally  arises  because  the 
values  of  dependent  variables  from  the  previous  process 
usually  do  not  satisfy  the  governing  equations  of  the  sub¬ 
sequent  process.  However,  if  there  is  no  outside  disturbance 
that  changes  the  differential  variables  (conservative  quan¬ 
tities),  the  consistent  initial  values  at  the  beginning  of  the 
succeeding  process  can  be  determined  from  the  ending  state 
of  the  previous  process  by  assuming  the  continuity  of  the 
differential  variables  [4].  This  can  also  be  easily  handled  by 
DAEIS  [7]. 

The  next  issue  of  serial  simulation  is  discrete  event 
detection  during  the  integration  of  a  continuous  process. 
For  serial  simulation,  the  ending  of  each  continuous  process 
is  due  to  the  occurrence  of  one  of  two  kinds  of  discrete 
events:  either  explicit  or  implicit  [4],  as  shown  in  Fig.  3.  An 


explicit  event  is  specified  by  the  independent  time  variable 
only,  which  is  actually  the  normal  working  mode  of  DAEs 
solvers  (i.e.,  specifying  the  ending  time).  However,  an 
implicit  event  is  based  on  some  dependent  variables,  which 
can  only  be  determined  during  the  integration  process.  The 
detection  of  discrete  events  requires  special  handling  in 
solving  DAEs,  e.g.,  checking  discrete  event  equations  at 
the  end  of  each  integration  step.  A  sign  change  of  the 
residual  evaluation  of  a  discrete  event  equation  represents 
the  occurrence  of  a  discrete  event.  Normally,  the  event 
location  is  determined  by  finding  the  root  of  interpolation 
polynomials  [2,4,8],  i.e.,  for  the  dependent  variables,  poly¬ 
nomials  can  be  constructed  on  an  interval  \tk_  i,  tk\  based  on 
the  past  solutions: 

Pi  =f(t,yf,yf~\yf~2,  ■■■)  (14) 

The  polynomial  equations  are  then  substituted  into  the  event 
definition  equations  (e.g.,  Eq.  (8)).  The  earliest  event  can 
then  be  located  by  combining  the  bisection  method  (deter¬ 
mining  the  time  interval  of  the  event)  and  Newton’s  algo¬ 
rithm  (finding  the  accurate  time  location  of  the  event)  [2]. 
With  the  determined  time  location  of  the  event,  the  values  of 
dependent  variables  at  that  moment  can  be  easily  obtained 
from  Eq.  (14). 

Few  DAE  solvers  are  capable  of  discrete  event  detection. 
DASRT,  an  extension  of  the  popular  DAEs  solver  DASSL 
[2],  is  used  in  this  work.  In  DASRT,  the  continuous  DAEs 
and  discrete  events  are  defined  in  separate  subroutines. 
When  an  event  is  detected  during  the  integration  process, 
DASRT  will  return  to  the  calling  program,  indicate  which 
event  occurred,  and  provide  the  state  of  dependent  variables 
when  the  event  happened. 

The  last  issue  of  serial  simulation  is  efficiency  and 
robustness  of  numerical  techniques.  The  modeling  of  elec- 


t0  tj  t2  t3  t4  t5 


Fig.  3.  Schematic  of  implicit  and  explicit  discrete  events. 
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trochemical  systems  usually  produces  hundreds  or  thou¬ 
sands  DAEs  (e.g.,  after  spatial  discretization  of  PDAEs). 
Solving  these  DAEs  normally  demands  significant  compu¬ 
tational  power.  In  addition,  an  electrochemical  system 
may  repeat  the  same  operations  for  many  times,  e.g.,  a 
secondary  battery  undergoes  thousands  of  cycles  of  charge/ 
open-circuit/discharge  processes.  It  can  take  days  or 
weeks  to  simulate  if  the  solving  efficiency  is  low,  which 
is  just  too  time-consuming  to  be  useful.  To  improve  the 
simulation  efficiency,  using  fast  hardware  is  a  simple  choice. 
However,  a  more  effective  approach  is  using  efficient 
numerical  techniques.  It  is  common  that  a  simulation 
taking  hours  to  finish  by  one  technique  (e.g.,  using  an  in- 
house  built  implicit  Euler  integration  code)  only  takes 
minutes  by  another  technique  (e.g.,  using  a  variable  order 
variable  step  size  solver  such  as  DASSL)  with  essentially 
the  same  results.  Besides  the  efficiency,  the  robustness  of 
the  numerical  techniques  is  another  critical  factor  in  the 
serial  simulation.  Coding  a  simulation  program  for  large 
DAEs  is  usually  a  nontrivial  task  and  it  is  even  harder  for 
serial  simulations.  Debugging  the  numerical  code  can  be 
terrible  if  the  robustness  of  the  numerical  techniques  utilized 
is  poor.  Actually,  simulation  failures  due  to  numerical 
algorithms  or  solvers  are  usually  impossible  to  detect  and 
correct.  Therefore,  a  simple  numerical  code  that  lacks 
sophisticated  error  controls  and  diagnoses  must  be  avoided, 
and  high-quality  professionally  built  numerical  packages 
should  be  used  instead.  There  are  many  efficient  and  robust 
solvers  for  DAEs,  some  of  them  are  available  free  from 
NETLIB  (http://www.netlib.org).  The  algorithm  of  the 
DAEs  solver  used  in  this  work,  DASRT,  is  based  on  the 
BDF  method  up  to  the  fifth  order.  At  each  integration  step, 
DASRT  uses  the  predictor  polynomial  to  predict  the  solution 
first  and  then  uses  the  corrector  polynomial  to  calculate  the 
final  results.  The  low-order  BDF  method  is  used  at  the  start 
of  integration  with  a  small  integration  step  size.  After 
enough  solution  points  have  been  obtained,  DASRT  will 
adjust  the  integration  step  size  and  BDF  order  based  on  the 
error  estimations.  There  are  many  details  involved,  which 
can  be  found  in  [2]. 

One  widely  used  ad  hoc  procedure  for  the  serial  simula¬ 
tion  is  to  use  a  fixed  time  step  integration  of  DAEs,  and  to 
check  if  some  conditions  (discrete  events)  are  triggered  after 
each  integration  step.  If  an  event  is  detected,  the  current 
process  is  terminated  and  the  final  state  of  the  process  is  used 
to  start  the  next  process.  However,  it  is  not  a  rigorous 
approach,  i.e.,  numerical  errors  due  to  the  inaccurate  termi¬ 
nation  of  the  previous  processes  will  accumulate  and  cause 
erroneous  simulation  results  for  the  later  processes.  Cer¬ 
tainly,  if  the  integration  time  step  is  chosen  to  be  sufficiently 
small,  e.g.,  10-5,  the  error  due  to  discrete  event  detection  of 
this  approach  will  be  insignificant,  but  the  simulation  effi¬ 
ciency  will  be  greatly  affected. 

With  the  consideration  of  above  issues,  a  procedure  for 
the  serial  simulation  of  electrochemical  processes  is  pro¬ 
vided  as  follows: 


1.  Construct  the  governing  equations  for  continuous 
processes.  If  PDAEs  are  obtained,  the  MOL  can  be 
used  to  transform  the  PDAEs  to  DAEs,  e.g.,  using  finite 
difference  approximations  for  spatial  derivatives.  For 
each  process,  identify  the  corresponding  discrete  events 
and  put  them  in  equation  form.  The  physical  constraints 
at  the  switching  between  processes  also  need  to  be 
determined,  e.g.,  which  conservation  law  will  apply. 

2.  Use  the  initialization  subroutine  DAEIS  to  provide  the 
consistent  initialization  for  the  DAEs  of  a  continuous 
process.  For  the  first  process,  initial  values  of  more 
accurately  known  variables  are  maintained  and  initial 
values  of  other  variables  are  determined  to  be  consistent 
with  those  variables.  When  switching  between  processes 
without  outside  disturbances,  values  of  differential 
variables  need  to  be  maintained  and  values  of  algebraic 
variables  are  determined  to  be  consistent  with  differ¬ 
ential  variables. 

3.  Apply  the  DAEs  solver  with  discrete  event  detection 
DASRT  to  solve  the  DAEs.  The  continuous  processes 
are  defined  in  one  subroutine  and  the  discrete  events  are 
specified  in  another  subroutine.  When  a  discrete  event  is 
detected  by  DASRT,  the  current  process  is  terminated. 
The  state  of  the  previous  process  at  the  discrete  event  is 
used  to  provide  the  starting  state  for  the  subsequent 
process. 

Although  the  above  procedure  seems  to  be  simple,  the 
complexity  of  the  relevant  numerical  algorithms  is  hidden 
inside  the  DAEIS  and  DASRT  solvers,  e.g.,  providing  a 
numerically  determined  Jacobian  matrix.  In  fact,  each  solver 
has  thousands  of  lines  of  carefully  constructed  code  to 
guarantee  the  efficiency  and  robustness  of  involved  algo¬ 
rithms.  To  use  the  above  procedure  effectively,  a  user  only 
needs  to  become  familiar  with  the  calling  protocols  of  the 
DAEIS  and  DASRT  solvers.  No  knowledge  of  the  algo¬ 
rithms  used  in  these  two  solvers  is  required.  Normally,  it 
takes  several  hundreds  lines  of  code  to  solve  a  medium  size 
model  in  the  combined  continuous/discrete  domain. 

For  the  given  example  problem,  less  than  200  lines  of 
fortran  code  were  used,  which  is  given  in  Appendix  A.  The 
DAEs  for  continuous  processes  are  given  in  the  subroutine 
DRES1.  The  discrete  events  are  defined  in  the  subroutine 
GR1.  The  residue  form  of  equations  are  used  in  both  DRES1 
and  GR1.  The  working  space,  solver  parameters,  and  equa¬ 
tion  definition  subroutines  are  shared  by  DAEIS  and 
DASRT.  The  main  subroutine  iterates  over  processes  and 
cycles.  It  is  obvious  from  the  code  that  the  details  of  solving 
algorithms  are  hidden  from  the  users  by  the  two  solvers, 
which  makes  the  usage  of  the  solving  procedure  much 
easier.  In  fact,  the  source  code  for  the  example  problem 
can  be  easily  adapted  to  handle  other  similar  simulations. 

The  results  of  the  simulation  are  shown  in  Figs.  4-8.  It  is 
obvious  from  Fig.  4  that  Processes  2  and  4  are  terminated  by 
the  explicit  discrete  events  of  time  specifications,  and  Pro¬ 
cesses  1  and  3  are  ended  by  the  implicit  discrete  events  of 
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t  (seconds) 


Fig.  4.  Simulated  potential  behavior  of  the  example  problem. 


cutoff  voltage  specifications.  Due  to  diffusion  limitations, 
there  are  nonuniform  concentration  distributions  of  NiOOH 
during  each  process  as  shown  in  Figs.  5-8.  Simulation  of  the 
cycling  like  that  shown  for  one  cycle  of  the  example 
processes  in  Fig.  4  can  be  done  for  thousands  of  cycles. 

From  the  results  of  the  example  problem,  it  is  clear  that 
much  more  information  can  be  obtained  from  the  serial 
simulation  of  multiple  processes  than  from  the  simulation  of 
a  single  process.  The  advantage  of  the  serial  simulation  is 
also  obvious:  investigations  that  are  impossible  to  conduct 
with  the  simulation  of  a  single  process  (e.g.,  the  open-circuit 
relaxation  behavior  after  a  charge  process)  can  be  easily 
achieved.  The  DAEs  resulting  from  the  modeling  of  elec¬ 
trochemical  systems  are  usually  more  complex  than  the 
given  example;  however,  the  same  numerical  solution  pro¬ 
cedure  applies.  In  fact,  the  procedure  has  been  utilized  in  the 
simulation  of  several  complex  models  of  battery  systems, 
including  a  nickel-hydrogen  cell  model,  a  nickel-metal 
hydride  cell  model,  and  a  lithium-ion  cell  model,  with 
satisfactory  results.  One  potential  application  of  serial  simu¬ 
lation  is  to  investigate  failure  mechanisms  by  simulating  the 
long  time  cycling  behavior  of  an  electrochemical  system 
with  a  model  including  degradation  processes. 


x 

Fig.  5.  Simulated  concentration  profile  of  the  example  problem  ( x  is  a 
dimensionless  coordinate)  for  Process  1  (charge  process). 
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x 

Fig.  6.  Simulated  concentration  profile  of  the  example  problem  ( x  is  a 
dimensionless  coordinate)  for  Process  2  (first  open-circuit  process). 


x 


Fig.  7.  Simulated  concentration  profile  of  the  example  problem  ( x  is  a 
dimensionless  coordinate)  for  Process  3  (discharge  process). 
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x 

Fig.  8.  Simulated  concentration  profile  of  the  example  problem  (x  is  a 
dimensionless  coordinate)  for  Process  4  (second  open-circuit  process). 

3.  Conclusions 

An  electrochemical  system  usually  operates  under  differ¬ 
ent  processes  consecutively.  To  simulate  the  corresponding 
behavior,  serial  simulation  of  many  processes  is  required, 
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which  is  demonstrated  with  the  simulation  of  charge/open-  ness  of  numerical  techniques.  It  is  shown  that,  with  two 

circuit/discharge  processes  of  a  nickel  hydroxide  electrode.  numerical  solvers,  DAEIS  and  DASRT,  serial  simulation  can 

Three  primary  numerical  issues  need  to  be  considered  in  be  easily  achieved.  Compared  to  the  simulation  of  a  single 

such  a  simulation:  consistent  initialization  of  a  process,  process,  the  serial  simulation  of  many  processes  allows  more 

accurate  termination  of  a  process,  and  efficiency  and  robust-  challenging  investigations  of  electrochemical  systems. 

Appendix  A.  Source  code  for  solving  the  sample  problem  with  DAEIS  and  DASRT 

C  Serial  simulation  of  a  thin  film  nickel  electrode. 

PROGRAM  NickelElectrode 

IMPLICIT  NONE 

INTEGER  NEQMAX,  MAXORD,  LRW,  LIW,  MAXPROC 
PARAMETER  (NEQMAX=501,  MAXORD=5,  MAXPROC=100) 

PARAMETER  (LRW  =  50+ (MAXORD+5) *NEQMAX+NEQMAX* *2 , 

&  LIW  =  20+NEQMAX+2*NEQMAX) 

DOUBLE  PRECISION  T,  TOUT,  RWORK (LRW) ,  RPAR(IOO) 

DOUBLE  PRECISION  Y (NEQMAX),  YPRIME (NEQMAX) ,  ATOL,  RTOL 
INTEGER  IWORK( LIW) ,  IPAR(2),  INFO(15),  IDID,  NEQ,  NG,  JROOT(2) 

EXTERNAL  DRES1,  DJAC1,  DDAEIS,  GRl 

INTEGER  I,  ICOUNT,  NCYCLE,  NPROC,  NODES,  NTOTCYCLE,  NTOTPROC 
DOUBLE  PRECISION  DTOUT ,  TEMPO,  THICKNI,  CNIOOHO,  XCURR (MAXPROC ) , 

&  XVCUT (MAXPROC) ,  XTCUT (MAXPROC) ,  RATENI,  RATE02 ,  DIFEH 

RTOL  =  l.D-5  !  Relative  Tolerance 

ATOL  =  l.D-9  !  Absolute  Tolerance 

OPEN (UNIT=10,STATUS=' unknown' , file=' op. txt' )  !  Input  File 

READ (10,*)  NODES  !  number  of  nodes 

READ(10,*)  DTOUT  !  time  interval  of  outputs  (seconds) 

READ(10,*)  TEMPO  !  ambient  temperature  (C) 

READ (10,*)  THICKNI  !  thickness  of  nickel  active  material  (cm) 

READ(10,*)  CNIOOHO  !  initial  NiOOH  concentration  (mole  fraction) 

READ (10,  * )  DIFFH  !  proton  diffusion  coefficient  (limitl  (cm2/s)) 

READ(10,*)  RATENI  !  exchange  current  density  of  nickel  reaction  (A/cm2) 

READ(10,*)  RATE02  !  exchange  current  density  of  oxygen  reaction  (A/cm2) 

READ(10,*)  NTOTCYCLE  !  total  number  of  cycles 

READ (10,*)  NTOTPROC  !  total  number  of  processes 

C  Read  operating  conditions  per  process 

DO  ICOUNT  =  1, NTOTPROC 

READ (10,*)  XCURR ( I COUNT ) , XVCUT ( ICOUNT ), XTCUT ( ICOUNT ) 

END  DO 
CLOSE (10) 

TEMPO  =  293.15  +  TEMPO 
RPAR (51)  =  96487. DO  !  constant 

RPAR ( 52 )  =  8.3141  !  constant 

RPAR (53)  =  TEMPO 

RPAR ( 54 )  =  3 . 4D0/92 . 71D0  !  maximum  proton  concentration,  mol/cm3 
RPAR (55)  =  DIFFH 
RPAR (56)  =  RATENI 
RPAR (57)  =  RATE02 

RPAR (58)  =  0.420  !  equilibrium  potential  of  nickel  reaction,  V 

RPAR (59)  =  0.303  !  equilibrium  potential  of  oxygen  reaction,  V 

RPAR (60)  =  THICKNI  !  thickness  of  nickel  active  material,  cm3 

NEQ  =  NODES  +  1 
I PAR ( 1 )  =  NODES 
NG  =  2 

OPEN(UNIT=12, STATUS=' unknown' ,FILE=' pot. txt' )  !  Output  File 

T  =  0.0D0  !  Initial  time. 

DO  NCYCLE  =  1, NTOTCYCLE 
DO  NPROC  =  1, NTOTPROC 

DO  I  =  2,  NEQ-2 

IWORK ( 2 0+NEQ+I )  =1  !  Y(I)  remains  initial  value 

IWORK (20+2*NEQ+I )  =  -1  !  YPRIME (I)  needs  to  be  determined 

END  DO 

IWORK ( 2  0  +NEQ+ 1 )  =  -1 


Y  ( I )  needs  to  be  determined 
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IWORK (20+2*NEQ+l )  =  0  !  YPRIME(I)  does  not  exist 

I WORK (20  +NEQ+NEQ- 1 )  =  -1 

IWORK (20+2*NEQ+NEQ-l)  =  0 

IWORK (20+NEQ+NEQ)  =  -1 

IWORK ( 20+2 *NEQ+NEQ)  =  0 

RPAR(l)  =  XCURR (NPROC) 

RPAR ( 2 )  =  T  +  XTCUT (NPROC) 

RPAR ( 3 )  =  XVCUT (NPROC) 

IF  (NCYCLE.EQ.l. AND. NPROC. EQ.l)  THEN 
DO  I  =  1,  NEQ-1 

Y ( I )  =  CNIOOHO  !  Initial  value. 

END  DO 

Y(NEQ)  =  0.40D0  !  Initial  value. 

END  IF 

DO  I  =  1,  NEQ 

YPRIME ( I )  =  -0.01  !  Initial  guess  values. 

END  DO 

DO  I  =  1,15 
INFO (I)  =  0 
END  DO 

CALL  DDAEIS (DRES1, NEQ, T, Y, YPRIME, INFO, RTOL, ATOL, IDID, 

&  RWORK, LRW, IWORK, LIW, RPAR, I PAR, DJAC1) 

150  CONTINUE 

TOUT  =  T  +  DTOUT 

CALL  DDASRT (DRES1, NEQ, T, Y, YPRIME, TOUT,  INFO, RTOL, ATOL, 

&  IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, D JAC 1 , GR1 , NG , JROOT) 

IF  (IDID  .LT.  0)  GO  TO  180 
IF  (IDID  .EQ.  4)  GO  TO  160 
CALL  RESULTS (NEQ, T,Y, YPRIME) 

GOTO  150 

160  CONTINUE 

IF  (IDID. EQ. 4 .AND. JROOT (1) .EQ.l)  THEN 
CALL  RESULTS (NEQ, T,Y, YPRIME) 

WRITE  (*,1000) 

1000  FORMAT (//IX, 1  Cut-off  time  reached',/) 

ELSE  IF  ( IDID . EQ. 4 . AND. JROOT (2 ) .EQ.l)  THEN 
CALL  RESULTS (NEQ, T,Y, YPRIME) 

WRITE  (*,1010) 

1010  FORMAT (//IX, '  Cut-off  voltage  reached',/) 

END  IF 
END  DO 
END  DO 

WRITE  (*,2000) 

2000  FORMAT (//IX, '  Simulation  Completed' , /) 

CLOSE (12) 

STOP 

180  CONTINUE 

WRITE  (*,2010) 

2010  FORMAT (//IX, '  Simulation  Failed ',/ ) 

CLOSE (12) 

STOP 

END 

SUBROUTINE  RESULTS (NEQ, T, Y, YPRIME) 

C  ================================================= 

IMPLICIT  NONE 
INTEGER  NEQ,  I 

DOUBLE  PRECISION  T,Y(*),  YPRIME (*) 

WRITE  (*,1020)  T, ( Y ( I ) , 1=1 , NEQ) 

WRITE  (12,1020)  T, (Y (I ) , 1=1, NEQ) 

1020  FORMAT (IX, 22 (E15.5) ) 

RETURN 

END 
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SUBROUTINE  DRES1 (T, Y, YPRIME, DELTA,  IRES,  RPAR,  IPAR) 

IMPLICIT  NONE 

INTEGER  IRES,  IPAR(*),  NEQ,  I 

DOUBLE  PRECISION  T,  Y ( * ) ,  YPRIME (*),  DELTA (*),  RPAR ( * ) 

DOUBLE  PRECISION  FJ1,  FJ2 ,  DELTAH 

DELTAH  =  RPAR (60) /FLOAT (IPAR(l) -1) 

NEQ  =  I PAR ( 1 )  +  1 

FJ1  =  RPAR (56)  *  (  (2.d0*  (l-Y(NEQ-l) ) )*exp(0.5* 

&  RPAR ( 51 ) /RPAR ( 52 ) /RPAR ( 53 ) * ( Y (NEQ) -RPAR ( 58 ) ) )  -  (2 . dO*Y (NEQ-1 ) ) 
&  *exp ( -0 . 5* RPAR ( 51 ) /RPAR ( 52 ) /RPAR (53) * ( Y (NEQ) -RPAR (58)  )  )  ) 

FJ2  =  RPAR (57)  * 

&  ( exp ( 1 . 0  *  RPAR (51) /RPAR ( 52 ) /RPAR(53) * (Y (NEQ) -RPAR(59)  ) ) 

&  -exp ( -1 . 0*RPAR ( 51 ) /RPAR(52) /RPAR (53) * ( Y (NEQ) -RPAR(59)  )  )  ) 

DELTA ( 1 )  =  (-3*Y (1) +4*Y(2)-Y (3) ) *RPAR ( 54 ) / DELTAH /2* RPAR ( 55 ) 

DO  I  =  2 , NEQ-2 

DELTA (I)  =  YPRIME (I) 

&  -  (Y(I+1) +Y(I-1) -2*Y(I) ) /DELTAH/DELTAH*RPAR(55) 

END  DO 

DELTA (NEQ-1 )  =  ( 3* Y (NEQ-1 ) -4 * Y (NEQ-2 ) +Y (NEQ-3 ) ) *RPAR ( 54 ) 

&  / DELTAH / 2  *  RPAR ( 55 )  -  FJ1/RPAR(51) 

DELTA  (NEQ)  =  FJ1  +  FJ2  -  RPAR(l) 

RETURN 

END 

SUBROUTINE  GR1  (NEQ,  T,  Y,  NG,  GROOT,  RPAR,  IPAR) 

IMPLICIT  NONE 

INTEGER  NEQ,  NG,  IPAR(*) 

DOUBLE  PRECISION  T,  Y(*),  GROOT (*)  ,  RPAR ( * ) 

GROOT (1)  =  T  -  RPAR (2) 

GROOT (2)  =  Y (NEQ)  -  RPAR (3) 

RETURN 

END 

SUBROUTINE  DJAC1 (T, Y, YPRIME, PD, CJ, RPAR,  IPAR) 

IMPLICIT  NONE 
INTEGER  I PAR ( * ) 

DOUBLE  PRECISION  T,  Y (*), YPRIME (*), PD (*), CJ, RPAR (*  ) 

C  DUMMY 

RETURN 
END 


The  input  file  ‘op.txf  for  the  program  is  given  below: 


20 

100.0 

25.0 

l.d-4 

0.05 

5 . d-12 

l.d-4 

l.d-10 

1 

4 

l.D-4 

0.0 

-1 .D-4 
0.0 

current ( 


number  of  nodes  in  the  nickel  active  layer 
time  interval  of  outputs  (seconds) 
ambient  temperature  (C) 

thickness  of  nickel  active  materials  (cm) 
initial  NiOOH  concentration  (mole  fraction) 
proton  diffusion  coefficient  (cm2/s) 
exchange  current  density  of  nickel  reaction  (A/cm2) 
exchange  current  density  of  oxygen  reaction  (A/cm2) 
number  of  repeating  times 
number  of  processes 

0.60  3600.0 

0.30  1800.0 

0.20  3600.0 

0.30  1800.0 

)  cut-off  Volt (V)  time  length  (sec) 
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